Applied Mound Detection

Chapter IV explains the development of workflow of the Master’s thesis and also the workflow itself. First the case study area (4.1) is presented with emphasis on its geomorphological characteristics, then a general workflow of automated analysis methods (4.2) is discussed. Following this, data pre-processing steps and methods are examined and the chosen method highlighted and the choice explained (4.3), similarly with the data preparation (4.4). Subsequently the data analysis methods, that is workflows used as paragons are discussed and the workflow of the Master’s thesis is explained (4.5).

The case study area

The case study area used in this Master’s thesis is a 180 km2 area around Marburg, in Marburg-Biedenkopf/Hesse, Germany (Figure 18). The research area lies in the central hessian region, mainly comprised of a densely populated settlement area (Marburg and its outskirts) in the center, to the West the Gladenbacher Bergland and to the East the mountainous ridge called Lahnberge (extending from north to south between Cölbe and Hassenhausen). West of the Lahnberge, which are made of triassic medium variegated sandstone extends the Amöneburg Basin with its loess landscape (Dobiat et al 1994, 7).

Location of the case study area of this Master's thesis in county Marburg-Biedenkopf/Hesse, Germany. Scale 1:100 000

Location of the case study area of this Master’s thesis in county Marburg-Biedenkopf/Hesse, Germany. Scale 1:100 000

The case study area was archaeologically investigated between 1983 and 1989 within the framework of a research project called “Urnenfelderzeitliche Forschungen im Marburger Raum” whose results were published in 1994 (Dobiat et al. 1994). This area constitutes the northernmost distribution of the South-German Urnfield Culture. Earlier studies led by the University of Marburg from the 1930’s outlined the Marburg Group (mainly coined by Müller-Karpe 1949) which based on the investigation between 1983 and 1989 can be broken down into multiple burial mound groups which extend over multiple hundreds of meters. Dobiat et al. 1994, 9 postulated that originally more than 250 burial mounds must have existed which formed groups of up to a dozen mounds and are mainly situated along ridges looking over the Amöneburg Basin. East of Marburg, North to South sprawl the three burial mound groups Botanischer Garten, Lichter Küppel and Stempel (Figure 19). According to the estimation of Dobiat et al. 1994, 31, the burial mound group of Botanischer Garten should have consisted of about at least 40 burial mounds while the burial groups at Lichter Küppel and Stempel of around 30 burial mounds each.

The burial mound groups ‘Botanischer Garten’, ‘Lichter Küppel’ and ‘Stempel’. Dobiat et al. 1994, Supplement 1.

The burial mound groups ‘Botanischer Garten’, ‘Lichter Küppel’ and ‘Stempel’. Dobiat et al. 1994, Supplement 1.

The extent of the case study area of this thesis was chosen intentionally based on the discussed area in Dobiat et al. 1994. After the 0.5 m DTM was calculated for the whole 180 km2 of the 20009/2010 LAZ files (see the workflow of the LiDAR data processing in 4.2), the attempt was made to identify all burial mounds in the Hillshade of the LiDAR derived 0.5 DTM by visual inspection in QGIS, based on the mapped burial mounds in Dobiat et al. 1994 (Map Supplement 1 & 2), also as a preparation to identify the Training Area and the Areas of Interests. Furthermore, apart from the already mentioned maps, Dobiat et al. 1994 also provides lists of burial mounds: Liste 1 contains the burial mounds of the three main burial groups East of Marburg (‘Botanischer Garten’, ‘Lichter Küppel’ and ‘Stempel’) and Liste 2 contains the burial mounds in the broader Marburg area (West to Marburg, Area of Interest 1). The exact coordinates can be found in Dobiat et al. 1994, 172-180. During the visual comparison of the provided lists of mapped burial mounds and the hillshade of the case study area in QGIS the following differences have been documented:

Comparison of the burial mound groups documented in Dobiat et al. 1994 and their visibility in the 2009/2010 LiDAR Hillshade. The sites which are marked with ? were not possible to identify in the Hillshade. NB: The burial mound groups 9 and 35 (the second mound) were only possible to identify using the Whitebox Multi-Scale Topographic Position Index derivative.
Site_ID Dobiat._et_al_1994 DTM05 AoI
1 13 mounds, some are still visible 12? visible 4
2 12 mounds 11-12? visible 4
3 5 + 2 mounds, not visible any more 5 + 1 + 1 visible 4
4 5 mounds 2 mounds visible? 4
5 8 mounds, some are still visible 9 visible 4/training area
6 5 mounds 3? visible 4
7 15 mounds, some are still visible 9 visible 4/training area
8 7 mounds 4 visible 4
9 2 mounds 1 visible in WB_MTPI 4/training area
10 8 mounds 4 visible? 4
11 2 mounds 2 visible 4
12 orig. 20 mounds, only few visible ~ 12 visible 4
13 1 mound 2 visible? 4
14 19 mounds 19 visible 4/training area
15 2 + 8 mounds 2? visible 4
16 4 flat mounds, hardly visible ? 3
17 13 mounds 10? visible 3
18 3 mounds destroyed 4
19 5 mounds destroyed? 3
20 ? mounds ? 3
21 7 mounds 3 visible 3
22 ~ 30 mounds? not clear if mounds ~ 20 visible 3
23 17 mounds ~ 11 visible 3
24 1 mound 2? mounds 3
25 ? mounds ? 3
26 6 mounds ? 3
27 3 mounds, destroyed? destroyed by road? 3
28 34 mounds, not all existent ~ 17? 3
29 at least 3 mounds 5? 3
30 1 mound destroyed by road? 2
31 ? mounds ? 2
32 1 mound? 1 visible 2
33 1 mound ? 2
34 1 mound 1+1 bit N ? 2
35 2 mounds in the field, ploughed in WB_MTPI visible 4/training area
49 2x2 mounds 2x2 mounds 1
51 1 + 2 mounds 1 + 2 mounds 1
61 2 mounds 2 mounds 1

The concept of the workflow for the automated analysis of this data set was developed with keeping in mind that the burial mounds marked by ‘? ‘ will not be detected with a high probability, if they are not traceable in the Hillshade, even when using profile tools in QGIS. It was clear that using a 180 km2 area to develop workflows on is overpowering the computational capacity of a regular Laptop (Tuxedo with Ubuntu 18.04.5 LTS, i7, 15.4 GB RAM, no extra graphic card) and even Desktop PC (Windows 10, i7, 36 GB physical RAM & 40.5 GB virtual RAM, NVIDIA GTX 1080), the data set was split up based on the investigation of the Hillshade of the case study area. First, all tiles were processed with 0.1 m resolution, but it was soon realized that any part of the workflow demanded exponentially more time per tile, than a tile with 0.5 m resolution, so the whole workflow is based on 0.5 m resolution.

The workflow developed was debugged in two instances: the workflow was first developed on the Training DTM, and then adapted to the Training Area. Then it was tested on five different Areas of Interests, which were chosen based on where the tumuli documented in Dobiat et al. 1994 are located. As Training DTM tile 3dm_32482_5618_1 was chosen (white box in Figure 20, left, with *‘_xyzirnc_ground_05’* suffix), because the burial mounds (Site IDs 5, 35 and 9, Figure 21) showed the biggest variability in size and were best preserved. Although it has to be said that even the Training DTM had surprises in store: Site ID 9 presented a quite a challenge to be identified - it was only possible to do so using the Multi-Scale Topographic Index implemented in Whitebox and accessed through R. The Training Area (red box in Figure 20) consisting of five 1x1 tiles (3dm_32482_5616_1_he, 3dm_32482_5617_1_he, 3dm_32482_5618_1_he, 3dm_32482_5619_1_he, 3dm_32483_5616_1_he, all with *‘_xyzirnc_ground_05’* suffix). Affiliated to these are the mound groups of Site ID 7 to the north, and Site ID 14, to the south (Figure 20, right).
The five Areas of Interest (AoI) are spatially disconnected were addressed clockwise: Area of Interest 1: Site IDs 49, 51 & 61 in List 2, Dobiat et al. 1994, 177-178. Figure 21 & 22. Area of Interest 2: Site IDs 30 - 34 in List 1, Dobiat et al. 1994, 175. See Figure 21 & 22. Area of Interest 3: ‘Botanischer Garten’ and ‘Lichter Küppel’, Site IDs 16 - 29 in List 1, Dobiat et al. 1994, 174-175. Figure 21 & 22. Area of Interest 4: ‘Stempel’, Site IDs 1 - 15 & 35 in List 1, Dobiat et al. 1994, 172-173 & 175. Figure 21 & 22. Apart from the burial mounds mapped and discussed in Dobiat et al. 1994, a few other burial mounds discernible in the Hillshade of the case study area were also included in the test dataset, that is AoI 5. Area of Interest 5: 4 tiles west to Gisselberg, displaying 8(?) merovingian burial mounds in the west corner, but also other mound-like structures can be determined, actually throughout the whole case study area.

Left:  Map and spatial relation of the Training DTM (white), the Training Area (red) and the five Areas of Interests (Area of Interest 1 = grey, Area of Interest 2 = blue, Area of Interest 3 = light green, Area of Interest 4 = grass green, Area of Interest 5 = orange). Right:The Training DTM, Training Area and the Area of Interests with the mapped burial mound sites in Dobiat et al. 1994. Scale 1:100 000Left:  Map and spatial relation of the Training DTM (white), the Training Area (red) and the five Areas of Interests (Area of Interest 1 = grey, Area of Interest 2 = blue, Area of Interest 3 = light green, Area of Interest 4 = grass green, Area of Interest 5 = orange). Right:The Training DTM, Training Area and the Area of Interests with the mapped burial mound sites in Dobiat et al. 1994. Scale 1:100 000

Left: Map and spatial relation of the Training DTM (white), the Training Area (red) and the five Areas of Interests (Area of Interest 1 = grey, Area of Interest 2 = blue, Area of Interest 3 = light green, Area of Interest 4 = grass green, Area of Interest 5 = orange). Right:The Training DTM, Training Area and the Area of Interests with the mapped burial mound sites in Dobiat et al. 1994. Scale 1:100 000

To understand the nature and the morphological characteristics of the burial mounds in this area and also to understand if any change has happened since the LiDAR survey in 2009/2010, the Site IDs 5 and 35 were inspected in the terrain. At that time Site ID 9 was not identified in the Hillshade. It has to be noted that for the main part of the Master’s thesis the QGIS plugin ‘VoGIS Profile Tool’ was used and only in a late phase of the work was the ‘Profile Tool’ QGIS plugin discovered and used, which displays minimal changes in the terrain a lot better than the previous tool. On May 10th in the afternoon the author and two colleagues (Simon Seyfried, a fellow student and Bjön Schmidt, now working in the Landesamt für Denkmalpflege in Wiesbaden) surveyed the area of the Training DTM. To identify the mounds, a polygon layer created in QGIS containing all identified mounds in the Training DTM (Site IDs 5 and 35) based on the visual inspection of the 0.5 m Hillshade was exported to OruxMaps on an Android mobile phone to use as orientation and localisation.

Screenshot of Oruxmaps with the burial mounds with Site IDs 35 (to the West in the field), 5 (to the East in the forest), 7 (to the North, not investigated in the field). Only Site IDs 5 and 35 were investigated which can be found in the area of Training DTM.

Screenshot of Oruxmaps with the burial mounds with Site IDs 35 (to the West in the field), 5 (to the East in the forest), 7 (to the North, not investigated in the field). Only Site IDs 5 and 35 were investigated which can be found in the area of Training DTM.

The first burial mounds to survey were Site ID 35. It was postulated from Dobiat et al. 1994, that there might be two mounds present (Table 1). Figure 22 implies that it is only from the texture of the terrain that it is possible to guess any height difference between the two mounds. This is verified by the profile. It is clearly visible that Site ID 35-1 is at most preserved up to around 38 cm in height. One of the surveyors is depicted in Figure 23 as a scale. Site ID 35-2 is on a slope, already probably ploughed away, and/or the separating line (the 25 cm protuberance) in the field also destroyed a part of the second mound (Figure 22).

Left: Location of burial mound group Site ID 35-1 and 2 in the Training DTM. Right: Profile of burial mound group, Site ID 35. Scale 1:700Left: Location of burial mound group Site ID 35-1 and 2 in the Training DTM. Right: Profile of burial mound group, Site ID 35. Scale 1:700

Left: Location of burial mound group Site ID 35-1 and 2 in the Training DTM. Right: Profile of burial mound group, Site ID 35. Scale 1:700

Björn as scale on top of burial mound Site ID 35-1. The descending terrain can be only guessed.

Björn as scale on top of burial mound Site ID 35-1. The descending terrain can be only guessed.

This was detected only when a Multi-Scale Topographic Position Index using the whitebox package was calculated (Figure 24). This derivative was not taken into account when building the workflow, because the whitebox R package (which is the R interface for the Whitebox GAT FOSS GIS software) did not work when working with the Tuxedo laptop and Windows 10 Desktop PC. Since then the bug was fixed in the next update. The Multi-Scale Topographic Position Index will be discussed more in depth later on. Generally it can be said that burial mounds which are situated in agricultural fields (that is in an intensively used area) are more likely to be detected by multispectral drone imagery and the calculation of vegetation indices than by LiDAR derived DTMs. Thus it was a pleasant surprise to see the traces of the two burial mounds of Site ID 35 (at least presumably also the second - and maybe a third?) in the Train DTM /Train Area/Area of Interest 4 using the Multi-Scale Topographic Position Index. The question is if burial mound Site ID 35-2 has eroded so much in the last 20 years (since 1994) that it is not possible any more to clearly detect it and only to try to guess it’s location. In Figure 24 the yellow areas delineate areas which are elevated in the micro-, meso- and macro scale, thus exaggerating the topography. On this basis also the area right to Site ID 35-2 might be arguable, but when we look at the extended profile it only with imagination to , and the section between 50 and 70 meters on the Y axis is most probable, but definitely eroded and disturbed by the field separator.

Left: The Multi-Scale Topographic Position Index calculated using the whitebox package in R, with the burial mound groups Site ID 35-1 and 35-2, possibly 35-3?. Right: Extended and intentionally positioned profile of burial mound group Site ID 35. Note the scale is different being an MTPI. Scale 1:700Left: The Multi-Scale Topographic Position Index calculated using the whitebox package in R, with the burial mound groups Site ID 35-1 and 35-2, possibly 35-3?. Right: Extended and intentionally positioned profile of burial mound group Site ID 35. Note the scale is different being an MTPI. Scale 1:700

Left: The Multi-Scale Topographic Position Index calculated using the whitebox package in R, with the burial mound groups Site ID 35-1 and 35-2, possibly 35-3?. Right: Extended and intentionally positioned profile of burial mound group Site ID 35. Note the scale is different being an MTPI. Scale 1:700

The burial mounds of Site ID 5(-1 to -9) are much more clearer to identify, as is visible from Figure 25. To make it easier to spot Site ID 5-9, the burial mound polygons are projected on the Hillshade with 25% transparency. Site ID 5 is the complete opposite of Site ID 35 - the mounds are extremely well discernable in the 0.5 m Hillshade, only to realize that when checking their profile, that their peak is mainly less than half a metre. Only Site ID 5-1 and 5 are outliers according to their height.

Left: Hillshade of Site ID 5. Right: Hillshade of Site ID 5 overlayed by transparent mound polygons. Scale 1:1500Left: Hillshade of Site ID 5. Right: Hillshade of Site ID 5 overlayed by transparent mound polygons. Scale 1:1500

Left: Hillshade of Site ID 5. Right: Hillshade of Site ID 5 overlayed by transparent mound polygons. Scale 1:1500

The approximative height measurements of the burial mounds are as follows: Site ID 5-1 ~ 120 cm Site ID 5-2 ~ 40 cm Site ID 5-3 ~ 38 cm Site ID 5-4 ~ 42 cm Site ID 5-5 ~ 50 cm Site ID 5-6 ~ 60 cm Site ID 5-7 ~ 50 cm Site ID 5-8 ~ 45 cm Site ID 5-9 ~ 40 cm

Below the profiles and the location of the profiles of the individual burial mounds are depicted for clarification (Figure 26-1 to 26-10.)

Left row: Hillshade of Site ID 5-1 to 9. Right row: Profile of Site ID 5 5-1 to 9 overlayed by transparent mound polygons. Scale 1:1500Left row: Hillshade of Site ID 5-1 to 9. Right row: Profile of Site ID 5 5-1 to 9 overlayed by transparent mound polygons. Scale 1:1500Left row: Hillshade of Site ID 5-1 to 9. Right row: Profile of Site ID 5 5-1 to 9 overlayed by transparent mound polygons. Scale 1:1500Left row: Hillshade of Site ID 5-1 to 9. Right row: Profile of Site ID 5 5-1 to 9 overlayed by transparent mound polygons. Scale 1:1500Left row: Hillshade of Site ID 5-1 to 9. Right row: Profile of Site ID 5 5-1 to 9 overlayed by transparent mound polygons. Scale 1:1500Left row: Hillshade of Site ID 5-1 to 9. Right row: Profile of Site ID 5 5-1 to 9 overlayed by transparent mound polygons. Scale 1:1500Left row: Hillshade of Site ID 5-1 to 9. Right row: Profile of Site ID 5 5-1 to 9 overlayed by transparent mound polygons. Scale 1:1500Left row: Hillshade of Site ID 5-1 to 9. Right row: Profile of Site ID 5 5-1 to 9 overlayed by transparent mound polygons. Scale 1:1500Left row: Hillshade of Site ID 5-1 to 9. Right row: Profile of Site ID 5 5-1 to 9 overlayed by transparent mound polygons. Scale 1:1500Left row: Hillshade of Site ID 5-1 to 9. Right row: Profile of Site ID 5 5-1 to 9 overlayed by transparent mound polygons. Scale 1:1500

Left row: Hillshade of Site ID 5-1 to 9. Right row: Profile of Site ID 5 5-1 to 9 overlayed by transparent mound polygons. Scale 1:1500

To sum up it can be said that most of the burial mounds of Site ID 5 were easily identified in the terrain, apart from Site ID 5-9. In Figures 26-9 and 26-10 traces of big agricultural machines are visible which makes sense, given that there is a concrete road in the vicinity as seen in Figure 27. The situation discernable in the Hillshade has apparently worsened in the last 10 years. Figure 27 illustrates the local topography of Site ID 5-9 (center illustrated with red stick) to demonstrate that if it weren’t for Oruxmaps and the shapefile of the burial mounds, the field identification of the mound would’ve been completely missed. This also suggests that other less well preserved or barely visible sites in the LiDAR data from 2009/10 might also be more eroded or even completely diminished. All this information has to be kept in mind when choosing automated analysis methods. Only because a certain method is working for burial mounds or mound-like archaeological objects in a certain region, it might not be true for other regions. Because topography in general, the micro-topography and especially geomorphology (and their way of erosion) play a very important role in how archaeological objects change with time and also how they are preserved in the present state when they are being investigated. In these cases a compromise has to be made to successfully detect barrows as will be demonstrated in the following sub-chapters.

Documentation of burial mound group Site ID 5-9 in the field.

Documentation of burial mound group Site ID 5-9 in the field.

The general workflow of Automated Analysis Methods

When undertaking automated analysis of remote sensing data, there are certain general steps to be followed - clearly outlining themselves in the surveyed studies and depending on what data type is at disposal. If DTM/DEM raster are available, Data Pre-Processing (step i) would be unnecessary. Also, some methods and workflows do not require data preparation (they directly use the DEM/DTMs) and start right away with Data Analysis (iii), which can be (as already detailed in Chapter 2) a combination of Geometric knowledge-based Analysis, GeOBIA and/or Machine Learning or single methods. In step iii the respective methods used in this thesis are investigated and discussed (iii-i and iii-ii). The last step (iv) is of course the verification of the results. A generalized workflow would look the following (see also Figure 28):

Generalized workflow of automated analysis of archaeological remote sensing datasets.

Generalized workflow of automated analysis of archaeological remote sensing datasets.

Data Pre-Processing

The LiDAR dataset from 2009/2010 of the case study area is processed by the state-of-the art point-cloud oriented lidR package (Roussel et al. 2020) which makes it easy to manipulate and analyse LiDAR datasets and is FOSS. Multiple functions, algorithms, methods and workflows presented in peer-reviewed literature are implemented and thus lidR enables users to make specialized workflows and thus functions as a toolbox (Roussel et al. 2020).

DTM generation by testing different algorithms

As a first step it was tested which ground point re-/classification and spatial interpolation should be used to generate a DTM which can capture all geomorphological details in the terrain that can be used to detect burial mounds. The processing of LiDAR data is a basic but very important part of the workflow. Still, it is only a means to an end so it will be kept brief (but in the necessary length to provide ideas to others who want to process LiDAR data for archaeological use using the lidR package) so that the focus can be oriented towards the detection of burial mounds.

Workflow of the LiDAR data processing.

Workflow of the LiDAR data processing.

For the comparison of the algorithms one random LAS tile was processed to understand the needed steps to generate at a DTM with as few artifacts and noise as possible, but detailed enough to still contain the information important for the Objects of Interest. A random tile was chosen, because when getting a point cloud data set we do not know yet where exactly our Objects of Interests are. The steps undertaken and tested are depicted and summarized in Figure 29. The complete code for this subchapter can be found in script 2_LiDAR_data_processing_tile_1.R. Here only relevant chunks are displayed.

4.3.1.0 Get to know the point cloud
LAS files store x,y,z for each point and many other information/attributes and this can claim a lot of memory from the PC.

names(LIDR_200910_1@data)
#[1] "X"                 "Y"                 "Z"                 "gpstime"
#[5] "Intensity"         "ReturnNumber"      "NumberOfReturns" #"ScanDirectionFlag"
#[9] "EdgeOfFlightline"  "Classification"    "Synthetic_flag"    "Keypoint_flag"
#[13] "Withheld_flag"     "ScanAngleRank"     "UserData"         "PointSourceID"

The ‘select’ argument enables you to choose between attributes/columns. The queries are: t - gpstime, a - scan angle, i - intensity, n - number of returns, r - return number, c - classification, s - synthetic flag, k - keypoint flag, w - withheld flag, o - overlap flag (format 6+), u - user data, p - point source #ID, e - edge of flight line flag, d - direction of scan flag

The ‘filter’ argument enables you to choose between the rows/points with specific attributes. The filter options can be accessed by: read.las(filter = "-help")

Note: when using the select and filter arguments with readLAS, it allows you to filter while reading the LAS file thus saving memory and computation time - it is the same as reading the LAS file and then filtering the point cloud but saving memory and computation time.

On this ground it was decided that the following arguments are needed: x,y,z, return number and number of returns, intensity, classification:

LIDR_200910_1_xyzirnc <- lidR::readLAS(lsLIDAR_2009_10[1], select = "xyzirnc")

4.3.1.1 Ground point re/classification The next step involves the re/classification of the ground points of the LiDAR point cloud. This serves as a foundation for the DTM generation. LAS point clouds already have a classification (as noted above) according to the ASPRS Society, which has explicit LAS specifications LAS specifications. Class 2 corresponds to the ground points. Two ground (re)classification algorithms are implemented in lidR: the Progressive Morphological Filter (PMF) and the Cloth Simulation Function (CSF). Both reclassify the classified ground point as unclassified before (re)classifying them. The aim of this chapter is to understand if the existing ground classification can be improved by reclassifying the point cloud and if the amount of testing is worth the time. Thus the PMF and CSF algorithms and the existing ground classification are compared if they can return an even more accurate model of the ground surface.

A) Ground Classification i LIDR_200910_1_ground

LIDR_200910_1_ground <- lidR::readLAS(lsLIDAR_2009_10[1], select = "xyzirnc", filter ="keep_class 2")

We can check the results by making a transect (2 x 500 m to be able do discern points) and plotting the cross section (green stands for class ground (2)):

Plot of the sections of the ground classification of the test tile.Plot of the sections of the ground classification of the test tile.

Plot of the sections of the ground classification of the test tile.

Cross sections were plotted for all point re/classifications, but the sections are best observed by zooming in. Thus the cross section plots are depicted in the Supplements (Chapter 9).

B) Progressive Morphological Filter (PMF) PMF (Zhang et al. 2003) generates a raster surface from the point cloud which then is processed by multiple morphological operations until stability is reached (Roussel et al. 2020, 4., [lidRbook] (https://jean-romain.github.io/lidRbook/gnd.html#pmf))

The workflow of the PMF algorithm.

The workflow of the PMF algorithm.

Multiple algorithms settings were tested. ws controls the window size of the moving window of the algorithm in meters. th is the threshold for point heights above the ground surface to be considered a ground return, also in meters.

iia LIDR_200910_1_xyzirnc_pmf

LIDR_200910_1_xyzirnc_pmf <- lidR::classify_ground(LIDR_200910_1_xyzirnc, algorithm = pmf(ws = 5, th = 3))

iib LIDR_200910_1_xyzirnc_pmf_th1 The height threshold was set to 1 m to avoid vegetation and outliers but include scattered ground points:

LIDR_200910_1_xyzirnc_pmf_th1 <- lidR::classify_ground(LIDR_200910_1_xyzirnc, algorithm = pmf(ws = 5, th = 1))

iic LIDR_200910_1_xyzirnc_pmf_th05 The height threshold was set to 0.5 m to avoid vegetation and outliers and to approach the ground surface:

LIDR_200910_1_xyzirnc_pmf_th05 <- lidR::classify_ground(LIDR_200910_1_xyzirnc, algorithm = pmf(ws = 5, th = 0.5))

iid LIDR_200910_1_xyzirnc_pmf_th005 The height threshold was set to 0.05 m to avoid vegetation and outliers and to get the ground surface as good as possible including possible backscatter of the LiDAR beam. When a height threshold of 0 was chosen, no ground points were (re)classified.

LIDR_200910_1_xyzirnc_pmf_th005 <- lidR::classify_ground(LIDR_200910_1_xyzirnc, algorithm = pmf(ws = 5, th = 0.05))

Comparing the section of iid with a threshold of 5 cm (0.05 m) with the sections of iic, with a threshold of 50 cm (0.5 m) is clearly visible that in the case of iid many of the actual ground points have not been classified as such because of the low threshold.

C) Cloth Simulation Function (CSF) CSF (Zhang et al. 2016) simulates a cloth which is dropped on the inverted point cloud and the ground points are classified by analyzing the interactions between the points of the cloth and the inverted surface (Roussel et al. 2020, 5; [lidRbook] (https://jean-romain.github.io/lidRbook/gnd.html#csf))

The workflow of the CSF algorithm.

The workflow of the CSF algorithm.

iiia LIDR_200910_1_xyzirnc_csf First the default settings were tested:

csf(sloop_smooth = FALSE,
     class_threshold = 0.5,
     cloth_resolution = 0.5,
     rigidness = 1L,
     iterations = 500L,
     time_step = 0.65)

Where sloop_smooth= TRUE reduces errors during post-processing; class_threshold is the distance to the simulated cloth to classify a point cloud into ground and non-ground; cloth_resolution is the distance between points in the cloth; rigidness stands for the rigidness of the cloth: 1 very soft (to fit rugged terrain), 2 medium, and 3 hard cloth (for flat terrain); time_step simulates the cloth under gravity. The default value is 0.65 and is suitable for most cases. (lidR vignette for the csf function).

LIDR_200910_1_xyzirnc_csf <- lidR::classify_ground(LIDR_200910_1_xyzirnc, algorithm = csf())

iiib LIDR_200910_1_xyzirnc_csf2 Because we are interested in depicting the ground surface in detail we are setting sloop_smooth= TRUE and leave the other parameter at default.

csf_1 <- csf(sloop_smooth = TRUE, class_threshold = 0.5, 
             cloth_resolution = 0.5, rigidness = 1, 
             iterations = 500, time_step = 0.65)
LIDR_200910_1_xyzirnc_csf2 <- lidR::classify_ground(LIDR_200910_1_xyzirnc, csf_1)

iic LIDR_200910_1_xyzirnc_csf3 In a last attempt rigidness was set to 2, because the region of Hessen is quite flat but still hilly and thus later we can compare which of the two settings yield better results - also compared to all of the ground segmentation methods.

4.3.1.2 Spatial Interpolation In the step before the LAS point cloud has been reclassified (by different methods). In this step the classified points are interpolated to create a terrain model.

Based on the specificities of the Objects of Interest in question, that is the barrows, it is clear that a resolution should be chosen which is detailed enough to also detect small mounds or destroyed ones by plowing, but which is not so high that it results in a big .tif file. Thus first resolutions between 0.5 and 0.05 meters have been tested (0.5, 0.2, 0.1 and 0.05). A DTM of 0.5 m resolution is between 12,000 and 13,000 KB. A DTM of 0.2 m resolution is between 73,000 and 75,000 KB and a DTM of 0.1 m resolution is between 270,000 and 280,000 KB. If calculating a DTM of 0.05 m = 5 cm resolution, the size of a tile of 1 km is more than 1 GB, which is way too much considering that the 2014 dataset consists of 209 1 km LAS tiles. First the resolution of 0.1 m was chosen, but it had to be realized that without HCP possibility, 0.5m is the best and fastest resolution and option a laptop and a quite good Desktop PC can go processing a rather large dataset.

Three interpolation algorithms are implemented in lidR: Triangulated Irregular Networks (TIN), Inverse Distance Weighing (IDW) and Kriging. In the following the tested settings are demonstrated only on the existing ground classification (LIDR_200910_1_ground).

A) Triangulated Irregular Networks (TIN)

This algorithm uses the Delaunay triangulation to generate a triangular irregular network from the classified point data. The resulting DTM presents an irregular surface based on non-overlapping triangles. The accuracy depends on the amount and density of points available in the point cloud. The more isolated the points are, the bigger the triangles get thus leading to abrupt elevation changes in the surface and thus to a distorted representation of the ground surface, resulting in edge artifacts. After testing the resolutions of 0.05, 0.1, 0.2 and 0.5 m, the default settings tin() were used with 0.5m resolution:

iva LIDR_200910_1_ground_tin05